



rm(list=ls())


library(coda)


## Load samples
load("model1_samples.Rdata")
load("model2_samples.Rdata")
load("model3_samples.Rdata")


## combine samples from 2 chain runs
m1.out <- runjags::combine.mcmc(list(m1.out[[1]],m1.out[[2]]))
m2.out <- runjags::combine.mcmc(list(m2.out[[1]],m2.out[[2]]))
m3.out <- runjags::combine.mcmc(list(m3.out[[1]],m3.out[[2]]))



## Posterior summaries
# note: differences between both DP spec. are minimal for betas
sum.coda(m1.out)
sum.coda(m2.out)
sum.coda(m3.out)


# Difference in estimates between N() and dt(4) RE priors
diff.dt <- m1.out[,5:20] - m2.out[,5:20]
sum.coda(diff.dt)

# Difference in estimates between N() and DP RE priors
diff.dp <- m1.out[,5:20] - m3.out[,6:21]
sum.coda(diff.dp)



# Plot differences for some selected variables

# Plot normal RE, T-RE, DP-RE
layout(matrix(c(5,5,1,2,3,4), ncol=2,byrow=TRUE), heights=c(.5, 2,2))
par(mgp=c(0, .35, 0))

par(mar=c(3.4, .8, .8 ,.8))

plot(1, xlim=c(0.78, 1.65), ylim=c(0,4.7),  type="n", axes=FALSE, xlab="", ylab="")
lines(density(m1.out[,"scale"], bw=0.02), col="gray30", lwd=1.3, lty=1)
lines(density(m1dt4.out[,"scale"], bw=0.02), col="gray30", lwd=2, lty=3)
lines(density(m1dp.out[,"scale"], bw=0.02), col="gray30", lwd=1.6, lty=2)
axis(1, tcl=-.3)
mtext(expression(lambda), 1, line=2, cex=1)
box()

plot(1, xlim=c(0.11, 0.34), ylim=c(0,15.2),  type="n", axes=FALSE, xlab="", ylab="")
lines(density(m1.out[,"gamma"], bw=0.02), col="gray30", lwd=1.3, lty=1)
lines(density(m1dt4.out[,"gamma"], bw=0.02), col="gray30", lwd=2, lty=3)
lines(density(m1dp.out[,"gamma"], bw=0.02), col="gray30", lwd=1.6, lty=2)
axis(1, tcl=-.3)
mtext(expression(phi), 1, line=2, cex=1)
box()

par(mar=c(3.4, .8, .8 ,.8))

plot(1, xlim=c(-.6, -.09), ylim=c(0,8), type="n", axes=FALSE, xlab="", ylab="")
lines(density(m1.out[,"beta[1]"], bw=0.02), col="gray30", lwd=1.3, lty=1)
lines(density(m1dt4.out[,"beta[1]"], bw=0.02), col="gray30", lwd=2, lty=3)
lines(density(m1dp.out[,"beta[1]"], bw=0.02), col="gray30", lwd=1.6, lty=2)
axis(1, tcl=-.3)
mtext(expression(beta[9]), 1, line=2, cex=1)
box()

plot(1, xlim=c(-.03, 0.54), ylim=c(0,7), type="n", axes=FALSE, xlab="", ylab="")
lines(density(m1.out[,"beta[11]"], bw=0.02), col="gray30", lwd=1.3, lty=1)
lines(density(m1dt4.out[,"beta[11]"], bw=0.02), col="gray30", lwd=2, lty=3)
lines(density(m1dp.out[,"beta[11]"], bw=0.02), col="gray30", lwd=1.6, lty=2)
axis(1, tcl=-.3)
mtext(expression(beta[10]), 1, line=2, cex=1)
box()

par(mar=c(.1, .1, .1 ,.1), mgp=c(0, .35, 0))
plot(1, xlim=c(-1,1), ylim=c(0,1), type="n", axes=FALSE, xlab="", ylab="")
legend(0, .4, c("Normal distributed", "t-distributed, df=4", "Mixture of dirichlet process"), col=c("gray30", "gray30", "gray30"), lty=c(1,3,2), lwd=c(1.6,2,1.6), bty="n", cex=1.1, xjust=.5, yjust=.45)




